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Impurity solvers play an essential role in the numerical investigation of strongly correlated elec- 
trons systems within the "dynamical mean field" approximation. Recently, a new class of continuous- 
time solvers has been developed, based on a diagrammatic expansion of the partition function in 
either the interactions or the impurity-bath hybridization. We investigate the performance of these 
two complimentary approaches and compare them to the well-established Hirsch-Fye method. The 
results show that the continuous-time methods, and in particular the version which expands in the 
hybridization, provide substantial gains in computational efficiency. 
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I. INTRODUCTION 

The numerical investigation of strongly correlated 
fermion systems on a lattice is a fundamental goal of 
modern condensed matter physics. Unfortunately, the 
standard Monte Carlo simulation technique is hampered 
by the sign problem P, , which leads to an exponen- 
tially growing statistical error with increasing system size 
or decreasing temperature. Exact diagonalization, on the 
other hand, is only possible for a small number of sites 
because of the exponential growth of the Hilbert space 
i- 

A useful and numerically feasible approach to treat 
fermionic systems in the thermodynamic limit is the so- 
called dynamical mean field theory. Development of this 
field started with the demonstration by Miiller-Hartmann 
and by Metzner and VoUhardt 0, Q that the diagram- 
matics of lattice models of interacting fermions simplifies 
dramatically in an appropriately chosen infinite dimen- 
sional (or infinite coordination) limit. This insight was 
developed by Georges, Kotliar and co-workers [IJtI who 
showed that if the momentum dependence of the elec- 
tronic self-energy may be neglected (E(p, w) 
as occurs in the infinite coordination number limit, then 
the solution of the lattice model may be obtained from 
the solution of a quantum impurity model plus a self- 
consistency condition. 

Quantum impurity models are amenable to numer- 
ical study. A well-established method is the Hirsch- 
Fye algorithm [s^, which employs a discrete Hubbard- 
Stratonovich transformation to decouple four-fermion 
terms. The method however requires the discretization 
of imaginary time into equal slices, which restricts it 
to relatively high temperatures. Further, analysis of 
the "Slater-Kanamori" interactions relevant to partially 
filled d-orbitals requires a very large number of Hubbard- 
Stratonovich fields Q . This further complicates the sim- 
ulations. 

Recently, a new class of continuous-time impurity 
solvers has been developed [l^, [ll|, . These diagram- 
matic QMC approaches rely on an expansion of the parti- 
tion function into diagrams and the resummation of dia- 



grams into determinants. A local update Monte Carlo 
procedure is then used to sample these determinants 
stochastically. Two complementary approaches have 
been formulated: the weak coupling method ^To| uses a 
perturbation expansion in the interaction part, while the 
hybridization expansion method [ll|, treats the local 
interactions exactly and expands in the impurity-bath 
hybridization. In the weak-coupling case, the determi- 
nantal formulation, which eliminates or at least greatly 
alleviates the sign problem, originates from Wick's theo- 
rem. In the hybridization expansion, when starting from 
a Hamiltonian formulation, the determinants emerge nat- 
urally from the trace over the bath states . 

Both continuous-time methods appear to provide con- 
siderable improvements over Hirsch-Fye. Our pur- 
pose here is to compare the performance of the two 
continuous-time solvers to each other and to Hirsch-Fye 
in an objective way. We use the algorithms and mea- 
surement procedures as proposed in Refs. 8, 10, 11], and 
focus on the accuracy with which physical quantities can 
be obtained in a DMFT calculation, for fixed CPU time 
in the impurity solver step. 

The rest of this paper is organized as follows: section 
II presents the needed formalism, section HI describes 
aspects of the measurement procedure, section IV gives 
results and section V is a conclusion. 

II. THEORY 
A. Model 

We concentrate in this work on the Hubbard model. 
For one band, the corresponding single-site impurity 
model is specified by the imaginary time effective action 

pP 

ScS = - drdr' V'cCT(r)_Fcr(T - r')4(T') 

-J dT^fi{n-^ +ni) -UuTni , (1) 

where n denotes the chemical potential and U the on- 
site repulsion. The hybridization function F describes 
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transitions into the bath and back and is related to the 
mean field function Qq by 0, [13] 



(2) 



The task of the impurity solver is to compute the Green 
function 



G{T-T') = -{T,c{Ty{T'))s^,,^ 



for a given hybridization function. 



TrTre-^-"c{T)cHT') 



TrTre-S-a 



(3) 



B. Hirsch-Fye impurity solver 

The algorithm of Hirsch and Fye Q requires a dis- 
cretization of imaginary time into N slices At = f3/N. In 
each time slice, the four-fermion term U n^ni is decoupled 
using a discrete Hubbard-Stratonovich transformation, 



-ArUir. 



^ s=±l 



(4) 



where the parameter A is defined as A = cosh(e'^'^'^/^). 
The Gaussian integral over the fermion fields may then 
be performed analytically, yielding an expression for the 
partition function of the form 



^det 



D, 



;-i_-p(si,...,SAf)£>g-i_j^(si,...,sjv) • (5) 



Here, Dg-i ^{si, sn) denotes the N x N matrix of the 
inverse propagator for a particular configuration of the 
auxiliary Ising spin variables Si, . . . , sjv [21 ■ The Monte 
Carlo sampling then proceeds by local updates in these 
spin configurations. Each successful update requires the 
calculation of the new determinants in Eq. ([5]) , at a com- 
putational cost of 0{N'^). 

The problem with this approach is the rapid (and, 
for metals, highly non-uniform) time-dependence of the 
Green functions at low temperature and strong interac- 
tions. The initial drop of the Green function is essen- 
tially - fr om which it follows that a fine grid 
spacing N ~ (3U is required for sufficient resolution. In 
the Hirsch-Fye community, N — j3U is apparently a com- 
mon choice, although we will see below that this number 
is too small and leads to significant systematic errors. As 
noted in Ref. [ll| a resolution of at least N = 5/3J7 is typi- 
cally needed to get systematic errors below the statistical 
errors of a reasonably accurate simulation. 



At half filling, the matrices Dg-i ^ and 



identical and it then follows immediately from Eq. ([5]) 
that the Hirsch-Fye algorithm under these conditions 
does not suffer from a sign problem. In fact, a closer 
analysis reveals that the sign problem is absent for any 
choice of /x [isj . 



C. Weak coupling expansion 

Recently, Rubtsov and co-workers proposed an entirely 
different approach for solving quantum impurity models 
[lo| . Their continuous-time method is a diagrammatic 
QMC algorithm which can be regarded as an extension 
of ideas originally introduced in Ref. [l^ to fermionic 
systems. The algorithm is based on a diagrammatic ex- 
pansion of the partition function in the interaction term 
and a stochastic sampling of the resulting diagrams (see 
Fig. [1]). More specifically, the action ([T]) is decomposed 
into a quadratic part 

5o - - / dTdT'y^C^{T)F^{T-T')cl{T') 

Jo 



-H I dT{ni +ni) 
Jo 



and an interaction part 



Sjj = U drn^ni. 



(6) 



(7) 



The weak coupling expansion of Z ~ TrTrC (So+Su) 
powers of U then reads 

Z = ^^-S^ I dn... dnTrTre~^° 
^ k\ J 

xn|(Ti)n|(Ti)...n|(Tfc)ri|(Tfc). (8) 

The trace over the fermionic degrees of freedom can now 
be performed analytically. Wick's theorem leads to 2kl 
terms whose combined weight is the determinant of the 
matrix product Dg^^i ^Ti, ...,Tk)Dgg^i{Ti, r^). The ele- 
ments of these k x k matrices are given by the mean field 
function defined in Eq. ^ 



DgoA^li ■■■,Tk){i,j) = Go,a{n ~ T.j). 

Finally, the partition function becomes 

k 



k\ 



dri . . . dTk det 



(9) 



(10) 



and the Monte Carlo sampling proceeds by local updates 
(random insertions/removals of vertices). At first sight, 
it appears that the term {—U)^ would lead to a bad 
sign problem for repulsive interactions. Rubtsov and co- 
workers found a way to get around this problem by re- 
defining the interaction term Sjj with a small positive 
constant a as 



Sf. ^ 'i 



dT 



{n^ir) +a){ni{T) - 1 - a) 
-(nt(T)-l-a)(nx(T)+a)l (11) 



and adjusting the quadratic term Sq in a way to compen- 
sate for this change [Tol|. However, this suppression of the 
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sign problem can lead to a probability distribution p„ (k) 
in the expansion order which features multiple maxima 
and the sampling of the different orders then requires 
some flat- histogram (or similar) method (isj . 



i) 




ii) 



FIG. 1: Illustration of the diagrams generated by the 
continuous-time impurity solvers, i) Weak-coupling method: 
third order diagram consisting of three vertices (dia- 
monds) I7n|(r)n|(r) linked by lines representing the function 
Go,cT{ji ~ Tj). ii) Hybridization expansion method: here, the 
orders of the diagrams for up- and down-spins can be differ- 
ent. Each creation operator cI{ts) (empty dot) is connected 
to an annihilation operator Co-(re) (full dot) by a line rep- 
resenting the hybridization function Faije — Ts). The black 
lines correspond to a particle number f , empty spaces to par- 
ticle number 0, so the overlaps between the lines for up- and 
down-spins yield the potential energy. In both approaches, 
the diagrams corresponding to different connecting t/o or F 
lines are summed up into a determinant and these determi- 
nants are sampled by a Monte Carlo procedure. 



D. Hybridization expansion 



Expanding the partition function Z = TrT^e (Sf+Sl) 
powers of F„ then leads to 



Ccr(Tl)Fo-(Tl - fl)4(fl) . . . 



(14) 



The individual terms in this series can have positive or 
negative sign, but as shown in Ref. [ll|, it is possible to 
express the combined weight of the fco-! diagrams corre- 
sponding to a given collection {cl{fi), Ca{Ti)}i=i^,,,^k„ of 
creation and annihilation operators as the determinant 
of a matrix Dp^a-, whose entries are the _F-functions, 

Dp.airi, . . . ,rfc„;fi, . . .,fkj{i,j) = F„{n - r,). (15) 

It can be proven in analogy to Ref. [l3| that this dia- 
grammatic formulation does not suffer from a sign prob- 
lem for models with density-density interactions'l!?!. The 
partition function finally becomes 



Z = 



r0 ffi 

a k^ -^0 -^^r 



dTk 



xc,{rl)cl{fl)...c^irncUrn, (16) 

where or denotes an upper integral bound which "winds 
around" the circle of length f3. If the last segment winds 
around, the sign is —1 and otherwise -1-1, whereas 
ST^ compensates for any sign change produced by the 
time ordering operator. The trace finds an easy and in- 
tuitive interpretation in terms of configurations of seg- 
ments marking the times where a particle of spin a is 
present [ll|. In such a representation, the /i-part of 5^ 
is determined by the total length of the segments while 
the interaction is given by the total overlap between seg- 
ments of opposite spin (see Fig. [T|) . 



A complimentary continuous-time algorithm (also il- 
lustrated in Fig. [T|) is obtained by expanding in the hy- 
bridization functions , while treating the chemical po- 
tential and interaction terms exactly. This approach has 
been worked out in Refs. [HI, [l2|- For the hybridization 
expansion, one decomposes the effective action ([1]) into 
the (non-local in time) hybridization part 

Sf = - dTdT'y^C„{T)F,{T-T')cl{T') (12) 

and the local part 

Sl — ^^J^ I driri] + ni) + U I drn^ni. (13) 
Jo Jo 



III. MEASURING THE GREEN FUNCTION 

A. Overview 

The diagrams obtained from the expansion of the par- 
tition function contain vertices or operators which are 
connected to each other by "bare" Green functions Qo 
or hybridization functions F. In order to measure the 
Green function, one needs a configuration with two "un- 
connected" operators c^(t) and c(t'). This may either be 
achieved by inserting such a pair of operators into a given 
configuration, or removing one of the Qq- or i^-functions 
from an existing diagram. The former approach has been 
implemented in the weak-coupling algorithm, and the lat- 
ter in the hybridization expansion algorithm. 
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B. Weak Coupling Expansion 

In the case of the weak couphng expansion the 
diagrams for the Green function G„{Tp — Tq) — 
— {^TCa{Tp)c]y{Tq)) SiTB obtained similarly to the expan- 
sion of the partition function (fTO|) , from the matrix Z?^' ^ 
containing an additional row and column: 



D 



Go, I 



GQ,aiTi - Tq) 
Go,a{Tp - Tq) 



(17) 



The relative weight of the Green function is given by 
the determinant ratio det Z?^' / det Dg^ ^ , which is com- 
puted using 



det D'^l = detDg^ 



\Go{tp 



E 



GoiTp ^ n){Dg^l),,goiT, - Tq) .(18) 



Hence the formula for measuring the Green function be- 
comes 

G{Tp,Tq) = gQiT.p-Tq)-/^go{Tp-T,)M,jgo{Tj-Tq)\, 



(19) 

where M = Dg^ and angular brackets denote the Monte 
Carlo average. Fourier transforming this formula yields 
a measurement formula in Matsubara frequencies, 

G{tun) = ao(«c^n)-(/3^'ao(«w»)'^AZ^je^""'"'-"^-^). 

(20) 

Both measurements can be performed directly during 
an update of the partition function, thereby reducing the 
computational effort for measuring the Green function 
from OiNM"^) to 0(NM), where N is the number of 
time slices or Matsubara frequencies and M the average 
matrix size [To| . 

The Matsubara Green function is required for the com- 
putation of the self-energy and the Hilbert transform, 
so measuring the Green functions directly in frequency 
space allows one to avoid the Fourier transformation from 
imaginary time to Matsubara frequencies. In the weak 
coupling algorithm the measurement in Matsubara fre- 
quencies appears as a correction to the (known) bare 
Green function C/o(*^n) which is suppressed by a factor 
of 1^^^ . For high frequencies, the errors converge very 
quickly and it is therefore possible to measure the high 
frequency behavior in a short time, before focusing on 
lower Matsubara frequencies for the rest of the simula- 
tion. This reduces the computational effort significantly. 



C. Hybridization Expansion 

Measurements in the hybridization expansion ap- 
proach can be performed by removing a hybridization 



function F{Tf — rj) connecting a pair of creation and 
annihilation operators. The contribution to the Green 
function at r = rf — is then given by the ratio 

(—1)*+^ AetD^p/ del Dp, where denotes the the ma- 
trix Dp with the row i and column j removed. Since 



{~iy+'deiD% = AetDF{Dp^)j, 



(21) 



it follows that the Green function can be measured from 
the inverse hybridization matrix M = Dp^ as 



T, T,- — T„- 



(22) 



where A(t, r') = sgn(T')^(r - t' - e{-T')(3). 

This yields 0{M^) (correlated) estimates for the Green 
function in one step. For r near or (3, the Green func- 
tion converges rapidly, whereas the variance of these mea- 
surement values is relatively large for r k, (3/2. The 
number of imaginary time slices does not influence the 
performance of the algorithm, and thus can be chosen 
arbitrarily large. It is therefore possible to adjust the 
resolution after the simulation according to the slope of 
the Green function and the noise in the measured data. 

Similar to the weak coupling case, formula could 
be Fourier transformed and G{iujn) measured directly. 
However, this approach is not advantageous, because the 
computation of the exponentials is very expensive com- 
pared to the rest of the simulation, and the binned values 
can easily be computed on a fine grid. Since the mea- 
sured values are in any case not a small correction to a 
known function, the hybridization expansion algorithm 
has more difficulties obtaining accurate results for the 
high-frequency behavior. 

As noted in Ref. [llj , the semi-circular density of states 
is a special case, where the self-consistency loop can be 
performed without any Fourier-transformation. The hy- 
bridization expansion results presented in the following 
section, however, were always obtained in the "conven- 
tional" way, by Fourier-transforming the Green function 
G(t) and extracting the self-energy. 



D. Susceptibilities 

In the weak coupling algorithm four-point functions 

Gi^j,(Ti,T2,T3,T4) = {TrcUTl)Ca{T2)cl,{T3)c^'{T4)) Can 

be measured by inserting operator pairs in analogy to 
the Green function ^Tf} . For same spin ct = cr' we then 
have to compute the determinant ratio of a matrix with 
two added rows and two added columns. For opposite 
spins one finds the product of two determinant ratios of 
the form (fT^ . 

In the hybridization expansion approach the 
four point functions are obtained by removing 
two hybridization lines, which leads to the mea- 
surement of J2ijkii^<^)i3i^'^')ki for CT ^ cr' and 
J2ijkiK^'y)ij(M'y)ki - {M„)aiMa)kj] for equal spins. 



5 



150 1 . , . , ^ 





u/t 

FIG. 2: Scaling of the matrix size with inverse temperature 
and interaction strength. Upper panel: temperature depen- 
dence for U/t = 4. In the case of Hirsch-Fye, the resolution 
N = l3U has been chosen as a compromise between reasonable 
accuracy and acceptable speed, while the average matrix size 
is plotted for the continuous-time solvers. Lower panel: de- 
pendence on U/t for fixed /3t — 30. The solutions for U < 4.5 
are metallic, while those for U > 5.0 are insulating. The much 
smaller matrix size in the relevant region of strong interactions 
is the reason for the higher efficiency of the hybridization ex- 
pansion method. 



Spin and charge susceptibilities, or more generally 
the density-density correlators can be obtained in- 
dependently and accurately (and with negligible 
computational effort) from the segment configurations. 
These segments, shown in the lower panel of Fig. [U 
represent the occupation of the orbital. 



IV. RESULTS 

A. Matrix size 

For all three algorithms, the computational effort 
scales as the cube of the matrix size, which for the 
Hirsch-Fye solver is determined by the time discretiza- 
tion At = f3/N and in the case of the continuous-time 



solvers is determined by the perturbation order k, which 
is peaked roughly at the mean value determined by the 
probability distribution p{k). In Fig. [21 we plot these 
matrix sizes as a function of inverse temperature /? for 
fixed U/t = 4 and as a function oi U/t for fixed /3t — 30. 
All our simulation results are for a semi-circular density 
of states with band-width it. 

It is obvious from the upper panel of Fig. [2] that the 
matrix size in all three algorithms scales linearly with /3. 
The Hirsch-Fye data are for N = f3U, which is apparently 
a common choice, although Figs. [3] and [5] show that it 
leads to considerable systematic errors. Thus, the grid 
size should in fact be chosen much larger {N > bfiU). 

While the matrix size in the weak coupling approach 
is approximately proportional to U /t, as in Hirsch-Fye, 
the {/-dependence of the hybridization expansion algo- 
rithm is very different: a decrease in average matrix size 
with increasing U/t leads to much smaller matrices in 
the physically interesting region 4 < U/t < 6, where the 
Mott transition occurs. The results in Fig. [2] and the 
cubic dependence of the computational effort on matrix 
size essentially explain why the continuous-time solvers 
are much more powerful than Hirsch-Fye and why the 
hybridization expansion is best suited to study strongly 
correlated systems. 

There is of course a pref actor to the cubic scaling, 
which depends on the computational overhead of the dif- 
ferent algorithms and on the details of the implementa- 
tion. Bliimer 16] has demonstrated substantial optimiza- 
tions of the Hirsch-Fye code and has in particular shown 
that extrapolating results at non-zero time step At to the 
At = limit considerably improves the accuracy. Of the 
continuous time codes investigated here, only the weak 
coupling results have been optimized. We estimate that 
similar modifications in the code for the hybridization ex- 
pansion algorithm would provide a speed-up of at least 
a factor of 10. However, the results presented here in- 
dicate large enough difference between the methods that 
the effects of optimization can be ignored. 



B. Accuracy for constant CPU time 

The three quantum Monte Carlo algorithms considered 
in this study work in very different ways. Not only are 
the configuration spaces and hence the update procedures 
entirely different, but also the measurements of the Green 
functions and other observables. 

In order to study the performance of the different im- 
purity solvers, we therefore decided to measure the accu- 
racy to which physical quantities can be determined for 
fixed CPU time (in this study 7h on a single Opteron 
244 per iteration). This is the question which is rele- 
vant to people interested in implementing either of the 
methods and avoids the tricky (if not impossible) task of 
separating the different factors which contribute to the 
uncertainty in the measured results. Because the vari- 
ance of the observables measured in successive iterations 
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of the self-consistency loop turned out to be considerably 
larger than the statistical error bars in each step, we de- 
termined the mean values and error bars using 20 DMFT 
iterations starting from a converged solution. 

The Hirsch-Fye solver suffers in addition to these sta- 
tistical errors from systematic errors due to time dis- 
cretization. These systematic errors are typically quite 
substantial and much larger than the statistical errors. In 
order to extract meaningful results from Hirsch-Fye simu- 
lations it is essential to do a careful (and time-consuming) 
At analysis [l^. The continuous-time methods are 
obviously free from such systematic errors if a sufficient 
number of time- or frequency points is used in the mea- 
surement of the Green function. 



1. Kinetic and potential energy 



0.52 




0.002 , , 0.004 



FIG. 3; Upper panel: kinetic energy Ekin ~ 
2t^ Jg dTG{T)G{-T) obtained using the three QMC im- 
purity solvers for U/t = 4.0 and f3t = 10, 15, ... , 50. The 
Hirsch-Fye simulations for At — 1/U (as in Fig. [Sjl yield 
systematically higher energies. The inset shows results 
obtained with the continuous-time solvers for f3t = 35, 40, 45 
and 50. Lower panel: potential energy U{n-^n{) for the same 
interaction strength. 



The kinetic energy, 

E^,^ = 2t^ / dTG{T)Gi~T), (23) 

^0 

shown in Fig. [3l was obtained from the imaginary time 
Green function by numerical integration. To this end 
we Fourier transformed the imaginary time Green func- 
tion and summed the frequency components including 
the analytically known tails. This turns out to be more 
accurate than the direct evaluation of equation by 
trapezoidal or Simpson rule. It is also more accurate than 
the procedure proposed in Ref. (l7j for the temperature 
and interaction range studied. 

We computed results for fixed U/t = 4 and temper- 
atures pt = 10, 15, ... , 50. In this parameter range the 
solution is metallic and we expect E'kin/^ oc (T/t)'^ at low 
temperature. The dominant contribution to £^kin comes 
from imaginary time points close to r = 0, /?. The accu- 
racy of the kinetic energy therefore illustrates how well 
the steep initial drop of G(r) can be resolved. 

The results from the continuous-time solvers agree 
within error bars, but due to the larger matrix size, the 
weak coupling algorithm can perform fewer updates for 
fixed CPU time and therefore the error bars are substan- 
tially larger (see inset of Fig. 

The Hirsch-Fye results are strongly dependent on the 
number of time slices used. Because of the cubic scal- 
ing of the computational effort with the number of time 
slices, at most a few hundred time points can be taken 
into account. This number is not sufficient to resolve the 
steep drop of the Green function at low temperature, and 
therefore the kinetic energy converges to values which 
are systematically too high. Extrapolation (e.g. ref. Q, 
[iHl) can be used to obtain values for At = and reduce 
these errors. However, various simulations at different 
At have to be performed in order to obtain an accurate 
estimate. For the kinetic energy we performed this ex- 
trapolation for f3t = 15, 20, 25. The error for (3 t = 20 
at At = after extrapolation is 10 times larger than 
the one we could obtain for the weak coupling algorithm, 
which is again around ten times larger than the one for 
the hybridization algorithm. 

We emphasize that for this particular case all three 
methods are sufficiently accurate that physically mean- 
ingful conclusions can be drawn; the differences, however, 
have clear implications for the extension of the method 
to more demanding regimes. 

In the lower panel of Fig. [3] we show the potential 
energy U{n-^ni) for U/t = 4, computed with the two 
continuous-time methods. In the hybridization expan- 
sion algorithm, the double occupancy can be measured 
from the overlap of the up- and down-segments. In 
the weak-coupling case, we used the relation U/2{{n^ + 
a){ni — 1 — a) + (ri| — 1 — a)(ri| +a)} = (k) / P (where (fe) 
is the average perturbation order), and an extrapolation 
to a — * 0. Both results agree within error bars and the 
hybridization expansion approach again yields the more 
accurate results. 
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FIG. 4: Lowest Matsubara frequency value for G for U/t — 
4.0, using measurements in both imaginary time and fre- 
quency space in the weak coupling case. The upper panel 
shows the Green function and the lower panel the relative er- 
ror on the measurement. Unlike in the Hirsch Fye algorithm 
there are essentially no systematic errors in the continuous 
time algorithms. In the case of the hybridization expansion 
algorithm, results for measurements in r and lu are plotted. 
Both measurements yield a similar accuracy at low frequency. 
The hybridization expansion algorithm gives very accurate re- 
sults and the error bars show no dependence on (3. This indi- 
cates that in the measured temperature range, two competing 
effects essentially cancel: the efficiency of the matrix updates 
which decreases at lower temperatures and the efficiency of 
the measurement procedure l|22[). which yields better results 
for larger matrix sizes. 



2. Green function and self energy 

The high precision of the hybridization expansion re- 
sults for the kinetic energy indicate that this algorithm 
can accurately determine the shape of the Green function 
near r = and /3. We now turn to the lowest Matsubara 
frequency component of the Green function, which is de- 
termined by the overall shape. We plot in Fig. [4] G(«a;o) 
for different values of p. The upper panel shows the re- 
sults obtained for the different continuous-time solvers 
and measurement procedures. They all agree within er- 



FIG. 5: Self-energy QmT,{iuJo) / ujq as a function of /3 for 
U/t = 4.0. The Hirsch-Fye results exhibit large discretiza- 
tion errors, while the continuous-time methods agree within 
error bars. The hybridization expansion method is particu- 
larly suitable for measuring quantities which depend on low- 
frequency components, such as the quasi-particle weight. 



ror bars. In the lower panel we plot the values of the 
error-bars. In the case of the weak-coupling expansion, 
both the measurement in r and the measurement in u 
produce about the same accuracy, which deteriorates as 
the temperature is lowered, due to the increasing matrix 
size. The error-bars from the hybridization expansion 
solver are much smaller and in the measured tempera- 
ture range remain about constant. Because the matrices 
at these values of U and /3 are very small, and the num- 
ber of measurement points in Eq. (|22p depends on the 
matrix size, the increase in computer time for updating 
larger matrices is compensated by a more efficient mea- 
surement. 

For the self-energy. 



(24) 



the Matsubara Green functions have to be inverted and 
subtracted. This procedure amplifies the errors of the 
self-energy especially in the tail region where C/o(*^n) and 
G{iLUn) have similar values. Fig. [5] shows 3mS(?a;o)/wo 
for U/t — A and several values of (3. This quan- 
tity is related to the quasi-particle weight Z w 1/(1 — 
^m'E,{iuJo)/ujo). Again, the Hirsch-Fye results show large 
systematic errors due to the time discretization and can- 
not be carried to low temperatures. The results from the 
continuous-time solvers agree within error-bars, but the 
size of the error bars is very different. The hybridization 
expansion approach yields very accurate results for low 
Matsubara frequencies in general. 

The advantage of measuring in Matsubara frequencies 
as opposed to imaginary time in the weak coupling algo- 
rithm becomes apparent for large ujn- Only the difference 
of G to the bare Green function Qq has to be measured 
in this algorithm. These differences decrease with l/a;„ 
for increasing a;„ and the estimate from Eq. (|20p is ex- 
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FIG. 6: Upper panel: low frequency region of the self-energy 
E(ia;) for U/t = 4.0, /3i = 45. Noise in tfie fiiglrer frequencies 
is clearly visible for the values measured in r, while the val- 
ues measured in ui in the weak coupling algorithm converge 
smoothly to the high frequency tail. Lower panel: high fre- 
quency region of the self-energy E(iti;) for U/t = 4.0, pt — 45. 
Noise in the higher frequencies is clearly visible for the values 
measured in r, while the values measured in lo in the weak 
coupling algorithm converge smoothly to the high frequency 
tail, limij^o '^{iiUn) = (7^(1 — n)n/{iu!„). 



tremely accurate at high frequencies, so that the tail of 
the self energy can be computed accurately. The mea- 
surements in imaginary time however have to be binned 
and Fourier transformed. While the high frequency tail 
can be enforced using correct boundary conditions for 
the cubic splines, there is a region of frequencies which 
starts much below the Nyquist frequency, where this in- 
troduces considerable errors (Fig.[S]). For 10 < LUn/t < 40 
and 500 imaginary time slices the values of show 
large errors before converging to the high-frequency tail 
enforced by the Fourier transformation procedure. The 
upper panel of Fig. [7] shows the difference between the 
two measurement approaches more clearly. 

The hybridization expansion algorithm starts from the 
atomic limit and thus does not get the high-frequency 
tail automatically right. Both a measurement in r and lu 
leads to relatively large errors at high frequencies. This 



FIG. 7: Intermediate frequency region of the self-energy T.{iu}) 
for U/t — 4.0, pt — 45. Noise in the higher frequencies is 
clearly visible for the values measured in r, while the val- 
ues measured in u in the weak coupling algorithm converge 
smoothly to the high frequency tail. 



noise again sets in at frequencies much below the Nyquist 
frequency, as illustrated by the results for 500 and 1000 
bins in the lower panel of Fig. [T] This noise is the con- 
sequence of the statistical errors in the Green function 
and can hence be reduced by running the simulation for 
a longer time (see Fig. [B]). However, Fig. [7] also shows 
that even for the shorter runs, the data remain accu- 
rate up to sufficiently large cun that a smooth patching 
onto the analytically known high-frequency tail appears 
feasible. Furthermore, since the hybridization expansion 
results in this section have all been obtained without any 
patching or smoothing and nicely agree with those from 
the weak-coupling solver, it seems that this uncertainty 
in the high-frequency tail is not a serious issue. 



C. Away from half filling 

We have tested both continuous time algorithms away 
from half filling, in a region where the half-filled model at 
zero temperature has a gap {U /t = 6, pt = 10) and in a 
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FIG. 8: Matrix sizes away from half filling: the matrix size 
decreases for the weak coupling algorithm, while the one for 
the hybridization expansion algorithm increases as one dopes 
a Mott-insulating state. 

region without gap [U/t = 3,l3t = 10, U/t = 2, /3t = 20). 
A comparison of the Green functions and self-energies 
has shown that both algorithms produce the same re- 
sult within numerical precision and are much faster than 
Hirsch-Fye. Both continuous time algorithms have no 
sign problem away from half filling (^13j, ,10,] ) and again 
the time needed to obtain a given accuracy is mostly de- 
termined by the size of the matrix. In the case of the 
weak coupling algorithm it decreases continuously away 
from half filling, while in the case of the hybridization ex- 
pansion the perturbation order first increases with doping 
if the half-filled model has a gap, and then decreases (see 
Fig.[8|). For all regions of parameter space tested, the hy- 
bridization expansion approach yields the smaller matrix 
sizes and is therefore substantially faster. The matrix 
sizes become comparable only in the limit of filled or 
empty bands. 

For the hybridization expansion algorithm, we have 
also computed the matrix size for U/t = 6 and much 
lower temperatures pt = 100, 200 and 400. These results 
showed that the perturbation order for a given filling re- 
mains proportional to /3, so that the shape of the curve 
remains the same as shown for /3t — 10 in Fig. [51 In 
particular this means that the formation of the "Kondo 
resonance" (which contains the physics of coherent low 
energy quasi-particlcs) in the slightly doped system at 
low temperatures does not lead to any dramatic change 
in the perturbation order. 

V. CONCLUSIONS 

We compared the performance of three different im- 
purity solvers: the Hirsch-Fye auxiliary field approach 



which for the last 15 years has been widely used in DMFT 
calculations, and two recently developed continuous-time 
algorithms, which are based on a stochastic sampling of 
an expansion of the partition function in the interac- 
tion and hybridization, respectively. Both continuous- 
time methods were found to be much more efficient than 
Hirsch-Fye for all relevant values of temperature, interac- 
tion strength and filling. Because the time-discretization 
in Hirsch-Fye simulations furthermore introduces a sys- 
tematic error which can be substantial and is difficult to 
estimate without a careful analysis involving several runs 
for different At, it makes sense (in most cases) to replace 
the method by one of the continuous-time solvers. 

The hybridization expansion leads to much smaller ma- 
trix sizes at intermediate and strong couplings, than the 
weak-coupling expansion. Hence it allows access to much 
lower temperatures. Our analysis has shown that the hy- 
bridization expansion is particularly powerful at calculat- 
ing low-frequency components of the Green function and 
quantities which are sensitive to these components, such 
as the quasi-particle weight. The method has more dif- 
ficulties capturing the high-frequency behavior correctly. 
In the weak-coupling approach, the algorithm perturbs 
around the non-interacting solution, which has the cor- 
rect high-frequency tail and if the Green function is mea- 
sured directly in frequency space, the intermediate and 
high Matsubara frequencies can be determined very accu- 
rately. Its power has also been demonstrated with recent 
applications to the Kondo lattice, multi orbital models 
and CDMFT 

We do not expect the noise in the high frequency com- 
ponents to be a serious problem, because the noise ap- 
pears at high enough frequencies that a smooth patching 
onto the analytically known tail seems feasible. However, 
the application of the hybridization expansion algorithm 
to a model with arbitrary density of states will be an im- 
portant test. Efforts to use it in the simulation of vana- 
dium oxide are under way and results will be presented 
elsewhere. 
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